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ABSTRACT 

The wall-matching methodology of Wilcox is modified to include a solid-wall, thermal-conduction 
model. This coupled fluid-thermal-structure model is derived assuming that the wall thermal-structure 
behavior is locally one-dimensional and that structural deformations, due to thermally induced stresses, 
are not significant. The one-dimensional coupled fluid-thermal-structure model is derived such that the 
wall temperature is removed as an independent boundary condition variable. The one-dimensional 
coupled fluid-thermal-structure model is also derived for the general case of an arbitrary mixture of 
thermally prefect gases and a wall of arbitrary thickness and conductivity by using a compressible, 
streamwise-pressure-gradient-corrected, wall-matching function and Fourier’s law of heat conduction. The 
resulting model was implemented in the VULCAN CFD code as a new boundary condition type. VULCAN 
was then used to simulate a two-dimensional Mach 6 wind tunnel facility nozzle flow to 
demonstrate/validate the one-dimensional coupled fluid-thermal-structure model. The nozzle internal-wall 
surface temperature and heat transfer distributions computed using the one-dimensional coupled fluid- 
thermal-structure model are compared to wall temperature and heat transfer distributions from an iterative 
multi-dimensional analysis obtained by coupling the VULCAN CFD code and the MSC/NASTRAN-thermal 
code. The one-dimensional coupled fluid-thermal-structure model analysis is shown to be very robust and 
in excellent agreement with the multi-dimensional iteratively coupled analysis. It is also shown that the 
one-dimensional analysis can be used as an initial guess for the multi-dimensional iteratively coupled 
analysis. 


INTRODUCTION 

The Hyper-X (1) and Advanced Space Transportation Programs (ASTP) <2) are currently 
investigating several airbreathing engine designs. The rocket based combined cycle (RBCC) and turbine 
based combined cycle (TBCC) engines of the ASTP and the Hyper-X program engine all incorporate 
scramjets as a component. These scramjet flowpaths are required to operate during the high Mach 
number (and high total enthalpy) portion of the vehicle flight trajectory. The operation of a scramjet engine 
at high total enthalpy levels causes the engine flowpath to be exposed to very high static temperatures, 
thereby presenting an extremely challenging thermal-structure design environment. At high total 
enthalpies, the near wall gas temperatures in the inlet and isolator exceed 2000 Kelvin and the combustor 
walls, fuel injector surfaces, flame holding devices as well as the nozzle walls are routinely exposed to 
gas temperatures in excess of 3000 Kelvin. The survival of structures utilizing available materials that are 
exposed to such high temperatures requires designs that are actively cooled. Furthermore, overall 
propulsion system thermodynamic-cycle efficiency and limitations on vehicle volume fraction constrain the 
vehicle design so as to require the use of the scramjet fuel as the structure coolant. 


The use of actively cooled structures encourages the development of analysis tools that provide 
accurate prediction of the structure heat loads and cooling requirements. However, engineering design 
tools for the analysis of the thermal environment have been predominately composed of either lower 
order methods, that attempt to approximately couple the fluid and thermal analysis, or higher order 
methods, that exclude coupling effects and use worst case scenarios to bracket the expected heat loads. 
Recent work performed as a part of the development of the National Propulsion Simulation System 
(NPSS) by Suresh (3) et al demonstrated a high-fidelity steady-state, weakly-coupled fluid-thermal- 
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structural solution methodology for scramjet flowpaths. This methodology iteratively coupled a CFD flow 
solver, Overflow* 4 ’, to a commercially available thermal-structural solver, ANSYS* 5 ’, with the assumption 
that the pressure and thermally induced deformations of the structure were negligible. The simulation that 
Suresh et at performed was an axi-symmetric analysis of the inlet of a GTX (3) engine module operating at 
Mach 4. Their fluid analysis began with an initial fluid-structure interface (wall) static temperature of 
277.77 Kelvin and their coupled solution ultimately converged to a wall temperature distribution that 
varied between 290 Kelvin and 370 Kelvin. Suresh et al observed that the coupling process converged 
relatively quickly (< 10 coupling iterations). However, their initial wall temperature guess was relatively 
close to the final temperature, which raises the concern that the number of coupling iterations could be 
unacceptable, when the initial guess is poor. Since the coupling process requires multiple time 
consuming (and expensive) multi-dimensional fluid and thermal structural analyses, it is desirable to 
develop methodologies that provide a better initial wall temperature boundary condition. It is also 
desirable to develop a methodology that would reduce the envelope where the use of an iterative, multi- 
dimensional, coupled fluid-thermal-structure analysis would be required. 


An approach that provides a better wall-temperature initial condition, and in some circumstances 
obviates the need for an iterative multi-dimensional coupled analysis, is to embed a one-dimensional 
thermal-structure analysis within the wall boundary condition used in the CFD analysis tool. The CFD 
code VULCAN* 6 ’, developed at NASA Langley Research Center specifically to analyze scramjet inlet, 
isolator, combustor and nozzle flows was chosen as the model development and demonstration platform 
due to its extensive use in scramjet flow analysis. VULCAN employs several techniques to reduce the 
total number of grid points that are required. One such technique is the use of the wall-matching function 
or “wall function” method of Wilcox* 7 ’. This wall-matching function method assumes that the first cell 
center off the wall is in the log law layer and that the wall shear stress is obtained by solving a modified 
form of the law of the wall equation. Wilcox’s method also utilizes a one-dimensional form of the energy 
equation to solve for the wall heat flux or to specify the wall temperature. The one-dimensional form of the 
near-wall energy equation can be exploited to provide a method through which the wall thermal structure 
analysis may be introduced. 


RESULTS AND DISCUSSION 

The wall matching method originally introduced by Wilcox* 7 ’ and later revised* 8 ’ begins with a 
modified statement of the compressible law of the wall that includes streamwise pressure gradient effects, 


“i‘ = u ,\t°ge(u T y 1 /vJ+KB + c 1 <l)(u T y i /vJ]/K i (1) 

where the streamwise pressure gradient term is defined as, 

0 = V 

and the wall friction velocity is, 



where p w ,v w , {dp/dx \ and z w are the wall density, kinematic viscosity, pressure gradient, and shear 

stress magnitude, respectively. The constants have the following values: c ; = -1.13, K = 0.41 and B = 
5. In addition, because the model is to be implemented in VULCAN (a cell centered code), the variables 
in equation (1) with a subscript of 1 refer to values taken from the first cell center adjacent to the wall. 
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Therefore, is the distance from the wall to the first cell center adjacent to the wall and u\ is the 
effective velocity tangent to the wall, in the first cell center adjacent to the wall. Wilcox (7) presents u\ as, 




where Pr f ,q w , C p ^ and T w are the turbulent Prandtl number, wall heat flux, effective specific heat at 

constant pressure (evaluated at the wall temperature) and the wall temperature, respectively, and is 
the velocity tangent to the wall, in the first cell center adjacent to the wall. 


Wilcox <7) presents the one-dimensional energy equation for the wall energy balance (generalized 
here for a frozen mixture of thermally prefect gases) as, 


C T - C T - -Pr 


1 2 

H U , 

2 1 


where 7j is the static temperature in the first cell center adjacent to the wall and C is the effective 
specific heat at constant pressure (evaluated at the static temperature in the first cell center adjacent to 
the wall). Equation (7) is then solved to yield an expression for T w , (which reduces to the adiabatic wall 

temperature if q w is set to 0), 

T = 

w 

or to yield an expression for q w , (which yields the wall viscous energy flux for a specified value of T w ), 
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MODIFICATION OF THE WALL MATCHING ENERGY EQUATION 


Fourier’s one-dimensional law of heat conduction* 9 * in differential form is 


q k =-k 


dT_ 

dy 


( 10 ) 


and is written for a planar wall, in discreet form as, 

<l- K 


l 


(t -t) 

V w e / 


(ii) 


w 


where k w is the thermal conductivity of the wall material, l w is the wall thickness and T e is the wall 
external temperature as shown in Figure 1. 
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Figure 1. - One-dimensional conduction through a plane wall with forced convection. 


Recognizing that, at steady state, the conduction heat transfer through the wall, q k , balances 
with the convection heat transfer from the gas to the wall, q w , yields, 

= (-It =-j-( T , r - T ,.) (12) 

'u 

which is solved for T w to give, 

T -T - — a 

W e 1 (13) 

w 

Substitution of equation (13) into equation (8) results in an equation that no longer requires that the 
flowpath gas-side wall temperature be known. The resulting equation can then be simplified to give, 
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( 14 ) 



Equation (14) implicitly couples the wall conduction and gas-side wall convection processes and is used 
instead of equation (9) when imposing the wall-matching function viscous-energy-flux boundary condition. 
Note that the superscript c has been added to indicate that the energy flux is now coupled to the wall 
conduction. The coupled wall-matching function process is then imposed on the gas-side flow solution 
through the following steps: 


1) An initial guess is created for the wall friction velocity, u T . 

2) Equations (2), (4), (5) and (6) are used to construct the components of equation (1) and are 
iteratively solved, to obtain u T , using Newton’s method. 

3) Equation (3) is used to compute the magnitude of the wall shear stress, z w . 

4) Equation (14) is used to compute the wall viscous energy flux due to coupled conduction and 
convection, q c w . 

5) The wall shear stress magnitude, z w , is used in the viscous momentum flux routines to 
impose the wall momentum flux. 

6) The wall energy flux due to coupled conduction and convection, q c w , is used in the energy 

equation viscous flux routine to impose the wall energy flux. 

7) Equation (13) is then used to provide the wall temperature for flow solution post-processing 
and visualization purposes. 


PROOF OF CONCEPT TEST CASE 

In order to test the wall-matching function coupling procedure, a relatively simple geometry that 
creates large, spatially varying, wall heat transfer loads is required. This requirement is satisfied by 
solving the flow through a two-dimensional convergent-divergent nozzle operating with high values of 
inflow total enthalpy and cooled walls. This nozzle flow solution is computed using two different 
approaches and the wall temperature and heat transfer distributions are compared. The first approach 
uses the one-dimensional coupled wall-matching function procedure implemented in the VULCAN code. 
The second approach uses a multi-dimensional coupled analysis similar to the one described by Suresh 
et at. However, in this instance, the VULCAN code, using an uncoupled wall-matching function method 
with a specified wall temperature distribution boundary condition, is iteratively coupled with the 
MSC/NASTRAN-thermal <10) analysis code. 


The test case selected was a two-dimensional nozzle design under consideration as a candidate 
to replace the existing three-dimensional Langley Arc-Heated Scramjet Test Facility' ' nozzle. The nozzle 
was computed using a single block two-dimensional grid consisting of 417 points in the i-direction and 65 
points in the j-direction. As shown in Figure 2, only one half of the nozzle height was modeled. The inflow 
boundary, located at the i-minimum boundary of the grid, was modeled as a subsonic inflow boundary by 
specifying the nozzle stagnation conditions. The outflow boundary, located at the i-maximum boundary of 
the grid, was modeled as an extrapolation boundary. The j-minimum boundary was modeled as a 
symmetry plane and the no-slip wall, located at the j-maximum boundary of the grid, was modeled using 
either a one-dimensional coupled wall-matching function or a specified temperature distribution wall- 
matching function. 
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Two-Dimensional Grid (417 x 65 Points) 



Figure 2. - Two-dimensional computational grid. 


The inflow boundary condition was modeled as a mixture of four chemical species and the flow 
chemistry was assumed to be chemically frozen throughout the expansion process. The tunnel plenum 
conditions were used to specify the inflow pressure and temperature as well as the gas composition and 
are presented in Tables 1 and 2, respectively. 


Stagnation Pressure 

35.384 atmospheres 

Stagnation Temperature 

2159.5 Kelvin 


Table 1. - Nozzle plenum/inflow pressure and temperature. 


Chemical Species 

n 2 

o 2 

Ar 

NO 

Mass Fraction 

0.7450 

0.2188 

0.0124 

0.0238 


Table 2. - Nozzle plenum/inflow gas composition. 


The nozzle calculations were run using the two-equation k-co turbulence model of Wilcox (8) 
with the flow assumed to be fully turbulent. A typical nozzle Mach contour plot resulting from these 
analyses, as presented in Figure 3, shows a well-designed nozzle flow with uniform outflow Mach number 
and an exit boundary layer thickness that is approximately 32 % of the nozzle exit height. The wall 
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Figure 3. - Typical nozzle Mach number contours. 


thermal analysis began by simulating the nozzle flow using a constant wall temperature of 500 K. 
Figure 4 presents the computed wall heat transfer distribution. The computed heat transfer 



Figure 4. - Computed nozzle wall heat transfer (Q) distribution for the isothermal wall case. 

reaches a peak of 1.5x10 7 Watts per meter squared at the nozzle throat. This heat transfer distribution 
was converted into an input file usable by MSC/NASTRAN-thermal and then a multi-dimensional (in this 
case 2-D) thermal analysis was performed to obtain the temperature distribution on the nozzle fluid-side 
wall. For simplicity, the nozzle wall was modeled as a constant thickness wall of 0.25 inches (6.35 mm.), 
with a uniform thermal conductivity of 14.4 Watts/Meter-Kelvin (a value typical of mild steel), and a 
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constant external temperature of 300 Kelvin. Once the multi-dimensional heat transfer problem was 
solved, the fluid-side wall temperature distribution was output from MSC/NASTRAN-thermal and 
converted into a VULCAN profile file (12) . The resulting temperature distribution illustrated a difficulty with 
using a poor initial temperature guess in the flow solver. The computed wall temperature was found to be 
significantly greater than the flow stagnation temperature. Ultimately, under-relaxing the temperature 
distribution obtained from the wall thermal analysis code during the initial iterations of the multi- 
dimensional coupling process alleviated this undesirable behavior. However, this under-relaxation 
increased the number of iterations required to converge the coupled analysis. 


Alternatively, the nozzle flow was computed using the one-dimensional coupled wall-matching 



CM 
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Figure 5. - Computed nozzle wall heat transfer (Q) distribution for the one-dimensional coupled 

wall matching function model. 



Figure 6. - Computed nozzle wall temperature distribution for the one-dimensional coupled wall 

matching function model. 
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function model given by equation (14). Figures 5 and 6 present the computed wall heat transfer and wall 
temperature distributions, respectively. Figure 5 shows that the nozzle throat computed wall heat transfer 
is approximately one fifth of the level developed by the iso-thermal wall temperature case. The reason for 
this is determined through inspection of the computed wall temperature distribution. Figure 6 shows that 
the computed wall temperature distribution reaches a maximum value of approximately 1600 Kelvin, 
which is approximately 5 times the assumed wall temperature used to perform the iso-thermal wall case. 
Note that the wall temperature for the majority of the nozzle length was computed to be approximately 
314 Kelvin, which is physically consistent with the assumed external surface temperature of 300 Kelvin. 


The computed wall heat transfer distribution obtained from the one-dimensional coupled wall- 
matching function model was then converted into an input file usable by MSC/NASTRAN-thermal and an 
iterative multi-dimensional analysis was initiated to study the feasibility of using the coupled one- 
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Multi-dimensional iterative method wall temperature distribution evolution. 


dimensional method as an initial guess for the multi-dimensional iterative method. The evolution of the 
wall temperature distribution during the iterative process is presented in Figure 7 for the initial one- 
dimensional coupled wall-matching function model method as well as 2 subsequent multi-dimensional 
coupled method iterations. The computed temperature distribution after iteration 2 can be seen to be very 
similar to the distribution computed using the one-dimensional coupled wall matching function model 
(iteration 0). Examination of a magnified view of the nozzle inflow and throat region shows that the 
differences between the iterations were less than 20 Kelvin at the nozzle throat, on the order of 80 Kelvin 
at the nozzle inflow, and less 5 Kelvin for the remainder of the nozzle. The differences between the one- 
dimensional and multi-dimensional methods at the nozzle inflow and throat are attributable to multi- 
dimensional effects. At the nozzle inflow the additional surface area of the leading edge of the wall that is 
in contact with the 300 Kelvin external temperature pulls more energy from the wall. This effect is 
completely ignored by the one-dimensional method. A similar observation can be made of the throat 
region of the nozzle. The throat axial temperature gradients on the fluid-side wall (and thus in the wall 
material) are large and cause significant conduction to take place. Again this effect is missed due to the 
nature of the one-dimensional coupled wall- matching boundary condition. Flowever, the overall 
agreement between iteration 0 and iteration 2 is excellent over the majority of the nozzle wall, 
demonstrating the utility of the one-dimensional method as an initial guess. This observation is reinforced 
by the heat transfer distribution after iteration 2, presented in Figure 8. When Figure 8 is compared 
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Figure 8. - Multi-dimensional iterative method wall heat transfer (Q) distribution after iteration 2. 


with Figure 5 the coupled one-dimensional wall matching model wall heat transfer distribution is seen to 
also be in close agreement with the iteratively coupled model wall heat transfer. 


SUMMARY AND CONCLUSIONS 

The wall matching methodology of Wilcox was modified to include a one-dimensional solid-wall 
thermal-conduction model. This coupled fluid-thermal-structure model was derived assuming that the wall 
thermal-structural behavior is locally one-dimensional and that structural deformations due to thermally 
induced stresses were not significant. The resulting model was implemented in the CFD code VULCAN 
as a new boundary condition type such that the wall thickness and thermal conductivity are either held 
constant or varied over the boundary surface. VULCAN was used to simulate a two-dimensional Mach 6 
wind-tunnel facility-nozzle flow. The predicted nozzle internal-wall surface temperature and heat transfer 
distributions were compared to the wall temperature and heat transfer distributions computed using a 
multi-dimensional analysis performed by iteratively coupling the VULCAN CFD code with the 
MSC/NASTRAN-thermal code in a manner similar to the approach of Suresh, et al. The one-dimensional 
analysis was found to be very robust and in excellent agreement with the multi-dimensional iteratively 
coupled analysis. It was also found that the one-dimensional analysis provides a much better initial guess 
for the multi-dimensional iteratively coupled analysis than the traditional iso-thermal guess, thereby 
reducing the number of iterations required to converge the multi-dimensional analysis. 
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